home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / roots / secant.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-05-14  |  2.9 KB  |  118 lines
  1. /* roots/secant.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Reid Priedhorsky, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* secant.c -- secant root finding algorithm 
  21.  
  22.    The secant algorithm is a variant of the Newton algorithm with the
  23.    derivative term replaced by a numerical estimate from the last two
  24.    function evaluations.
  25.  
  26.    x[i+1] = x[i] - f(x[i]) / f'_est
  27.  
  28.    where f'_est = (f(x[i]) - f(x[i-1])) / (x[i] - x[i-1])
  29.  
  30.    The exact derivative is used for the initial value of f'_est.
  31.  
  32. */
  33.  
  34. #include <config.h>
  35.  
  36. #include <stddef.h>
  37. #include <stdlib.h>
  38. #include <stdio.h>
  39. #include <math.h>
  40. #include <float.h>
  41.  
  42. #include <gsl/gsl_math.h>
  43. #include <gsl/gsl_errno.h>
  44. #include <gsl/gsl_roots.h>
  45.  
  46. #include "roots.h"
  47.  
  48. typedef struct
  49.   {
  50.     double f;
  51.     double df;
  52.   }
  53. secant_state_t;
  54.  
  55. static int secant_init (void * vstate, gsl_function_fdf * fdf, double * root);
  56. static int secant_iterate (void * vstate, gsl_function_fdf * fdf, double * root);
  57.  
  58. static int
  59. secant_init (void * vstate, gsl_function_fdf * fdf, double * root)
  60. {
  61.   secant_state_t * state = (secant_state_t *) vstate;
  62.  
  63.   const double x = *root;
  64.  
  65.   GSL_FN_FDF_EVAL_F_DF (fdf, x, &(state->f), &(state->df));
  66.   
  67.   return GSL_SUCCESS;
  68.  
  69. }
  70.  
  71. static int
  72. secant_iterate (void * vstate, gsl_function_fdf * fdf, double * root)
  73. {
  74.   secant_state_t * state = (secant_state_t *) vstate;
  75.   
  76.   const double x = *root ;
  77.   const double f = state->f;
  78.   const double df = state->df;
  79.  
  80.   double x_new, f_new, df_new;
  81.  
  82.   if (state->df == 0.0)
  83.     {
  84.       GSL_ERROR("derivative is zero", GSL_EZERODIV);
  85.     }
  86.  
  87.   x_new = x - (f / df);
  88.  
  89.   f_new = GSL_FN_FDF_EVAL_F(fdf, x_new) ;
  90.   df_new = (f_new - f) / (x_new - x) ;
  91.  
  92.   *root = x_new ;
  93.  
  94.   state->f = f_new ;
  95.   state->df = df_new ;
  96.  
  97.   if (!finite (f_new))
  98.     {
  99.       GSL_ERROR ("function not continuous", GSL_EBADFUNC);
  100.     }
  101.  
  102.   if (!finite (df_new))
  103.     {
  104.       GSL_ERROR ("function not differentiable", GSL_EBADFUNC);
  105.     }
  106.       
  107.   return GSL_SUCCESS;
  108. }
  109.  
  110.  
  111. static const gsl_root_fdfsolver_type secant_type =
  112. {"secant",                /* name */
  113.  sizeof (secant_state_t),
  114.  &secant_init,
  115.  &secant_iterate};
  116.  
  117. const gsl_root_fdfsolver_type  * gsl_root_fdfsolver_secant = &secant_type;
  118.